Modelos lineales generalizados (GLMs)

Introducción

  1. Definición

  2. Estructura

    2.1. Función de vínculo

    2.2. Familia de distribución de errores

1. Definición

Supuestos modelos lineales:

  1. Los errores se distribuyen normalmente

  2. La varianza (residual) es constante

  3. La variable respuesta se relaciona linealmente con la(s) variable(s) independiente(s), si alguna de estas variables es continua

Uno o varios supuestos no se cumplen.

¿Qué hacer?

GLMs: extensión modelos lineales especificando estructuras no normales de errores y varianzas no constantes.

Modelos lineales caso particular de GLMs con distribución de errores normal y varianza constante.

Relación varianza y media (Cayuela y de la Cruz, 2022)

Tipos de variables que por definición violan supuestos modelos lineales:

  • conteo de casos (p.ej. abundancia de especies)

  • conteo de casos como proporciones (p.ej. porcentaje plántulas muertas en un experimento)

  • respuesta binaria (p.ej. usado o no usado, presente o ausente)

2. Estructura

3 componentes principales:

1. El predictor lineal

2. La función de vínculo

3. La estructura de los errores

Estructura modelos lineales:

yi = β0 + β1x1 + … + β1xn + Ɛi

Ɛi ~N(0,σ2)

Estructura GLMs:

f(yi) = β0 + β1x1 + … + β1xn + Ɛi

Var(y) ~ f(µ)

2.1. La función de vínculo f(y)

Transformación de la variable y buscando linealizar relación y con x y/o transformando en continuas variables que no lo son.

Por tanto, cambia la naturaleza de la variable respuesta.

Determinadas funciones de vínculo deben implementarse con determinados tipos de variables respuesta

Función de vínculo (Cayuela y de la Cruz, 2022)
2.2. La familia de distribución de errores

  1. Poisson: varianza aumenta linealmente con la media. Para respuestas tipo conteo.

  2. Binomial negativa: extensión de Poisson. Para variables conteo con sobredispersión.

  3. Binomial: para proporciones y datos presencia/ausencia.

  4. Gamma: para datos continuos con coeficiente de variación constante.

Ajuste GLM en R

Función glm()

Argumento family:

  • binominal(link = “logit”)
  • gaussian(link = “identity)
  • Gamma(link = “inverse”)
  • poisson(link = “log”)

Ejemplo con datos tipo conteo

Estudio del efecto de la ganadería en la regeneración de araucaria. Muestreo de número de plantas de araucaria en dos tipos de ganadería: pequeños propietarios, con uso más intensivo de la tierra, y grandes empresas forestales, que mantienen algunos remanente de bosque nativo. Actividad ganado medida por número de bostas.

Variables:

  • Número de plántulas de araucaria

  • Número de bostas (proxy actividad del ganado)

  • Dos tipos de ganadería: pequeños propietarios y grandes empresas

library(ADER)
data(ara)
str(ara)
'data.frame':   26 obs. of  3 variables:
 $ seedlings: int  7 10 0 5 54 0 30 14 1 2 ...
 $ dungs    : int  5 0 330 960 115 2330 475 10 2000 130 ...
 $ property : Factor w/ 2 levels "Campesino","Empresa Forestal": 2 2 2 2 2 2 2 2 2 1 ...

Hipótesis:

H0: β = 0 (pendiente entre plántulas y número de bostas es 0)

H0: µ1 = µ2 (número plántulas en pequeños propietarios y empresas sin diferencias)

H0: β1 = β2 (relación entre plántulas y bostas es igual en pequeños propietarios y empresas)

Exploración gráfica de los datos

library(ggplot2)
ggplot(ara, aes(x = dungs, y = seedlings)) +
  geom_point(color = "steelblue", size = 2) +
  facet_wrap(~ property) +
  labs(
    title = "Relación nº plántulas y nº bostas",
    x = "Nº bostas",
    y = "Nº plántulas"
  ) +
  theme_minimal()

Construcción modelo

glm_arau1 <- glm(seedlings ~ dungs*property, ara, family = "poisson")
anova(glm_arau1, test = "Chi")
Analysis of Deviance Table

Model: poisson, link: log

Response: seedlings

Terms added sequentially (first to last)

               Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                              25     347.56              
dungs           1   30.967        24     316.59 2.624e-08 ***
property        1   80.598        23     236.00 < 2.2e-16 ***
dungs:property  1   37.974        22     198.02 7.170e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Devianza explicada por el modelo

Medida de cantidad de variabilidad que explica el modelo (similar R2 en modelos lineales)

D2 = 1 - Dmodelo/Dnulo

print(D2 <- 1 - (198.02/347.56))
[1] 0.4302566
summary(glm_arau1)

Call:
glm(formula = seedlings ~ dungs * property, family = "poisson", 
    data = ara)

Coefficients:
                                Estimate Std. Error z value Pr(>|z|)    
(Intercept)                     2.497423   0.135408  18.444  < 2e-16 ***
dungs                          -0.010211   0.001789  -5.706 1.15e-08 ***
propertyEmpresa Forestal        0.576261   0.170467   3.380 0.000724 ***
dungs:propertyEmpresa Forestal  0.009041   0.001803   5.015 5.31e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 347.56  on 25  degrees of freedom
Residual deviance: 198.02  on 22  degrees of freedom
AIC: 282.3

Number of Fisher Scoring iterations: 5

Fórmula: log(y) = 2.497 + 0.576 Empresa + (-0.010 + 0.009 Empresa) Bostas

Supuestos del modelo

par(mfrow=c(2,2))
plot(glm_arau1)      
library(aods3)
gof(glm_arau1)      
D  = 198.022, df = 22, P(>D) = 2.772377e-30
X2 = 248.4275, df = 22, P(>X2) = 2.969764e-40

Modelo binomial negativo

library(MASS)
glm_arau2 <- glm.nb(seedlings ~ dungs*property, ara)
summary(glm_arau2)      

Call:
glm.nb(formula = seedlings ~ dungs * property, data = ara, init.theta = 1.13715532, 
    link = log)

Coefficients:
                                Estimate Std. Error z value Pr(>|z|)    
(Intercept)                     2.152593   0.362295   5.942 2.82e-09 ***
dungs                          -0.005663   0.002468  -2.294   0.0218 *  
propertyEmpresa Forestal        1.113371   0.557993   1.995   0.0460 *  
dungs:propertyEmpresa Forestal  0.003945   0.002528   1.561   0.1186    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(1.1372) family taken to be 1)

    Null deviance: 49.506  on 25  degrees of freedom
Residual deviance: 28.040  on 22  degrees of freedom
AIC: 155.22

Number of Fisher Scoring iterations: 1

              Theta:  1.137 
          Std. Err.:  0.384 

 2 x log-likelihood:  -145.218

Supuestos del modelo

par(mfrow=c(2,2))
plot(glm_arau2)      
gof(glm_arau2)      
D  = 28.0403, df = 22, P(>D) = 0.1743489
X2 = 24.0659, df = 22, P(>X2) = 0.343783

Interpretación del modelo

library(car)
Anova(glm_arau2, type = "III")  
Analysis of Deviance Table (Type III tests)

Response: seedlings
               LR Chisq Df Pr(>Chisq)   
dungs            7.9871  1   0.004711 **
property         4.3145  1   0.037789 * 
dungs:property   3.5106  1   0.060977 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Predicciones del modelo

library(visreg)
visreg(glm_arau2, "dungs", by ="property", overlay = T, scale = "response", ylab = "Nº plántulas", xlab = "Nº bostas")

Ejemplo con distribución de errores binomial

2 variantes de estos modelos:

  • Variable con valores binarios (1/0): modelos de regresión logística

  • Variables con valores porcentuales (nº éxitos / nº de ensayos)

Función de vínculo: logit

log(p/(1-p))

logit(y) = ln(p/(1-p)) = β0 + β1x1 + … + β1xn

Curva sigmoidea entre 0-1

Ejemplo: Estrategia reproductora especie de alga invasora

Estrategia reproductora de un alga invasora con dos especies nativas.

¿Probabilidad de encontrar estructuras reproductoras es diferente entre la especie invasora y las nativas?

Hipótesis: la especie invasora tendrá mayor probabilidad de reproducirse que las nativas

Variable respuesta (Estado):

 - 0: sin estructuras reproductoras
 - 1: con estructuras reproductoras

Variables predictoras:

 - talla del individuo (cm)
 - especie (Sp)
 
library(ADER)
data(algas)
str(algas)
'data.frame':   359 obs. of  3 variables:
 $ Sp    : Factor w/ 3 levels "Fra","Tom","Ver": 3 3 3 3 3 3 3 3 3 3 ...
 $ Long  : num  1 2 4 4 4 4.5 4.5 4.5 5 5 ...
 $ Estado: int  0 0 0 0 0 0 0 0 0 0 ...

Exploración de los datos

ggplot(algas, aes(x = Long, y = Estado, color = Sp)) +
  geom_point(alpha = 0.4, position = position_jitter(height = 0.05)) +
facet_wrap(~Sp) +
  labs(
    x = "Talla individuo (cm)",
    y = "Estructura reproductora",
  ) +
  theme_minimal()

Ajuste del modelo

glm_algas1 <- glm(Estado ~ Long*Sp, algas, family = "binomial")
anova(glm_algas1)
Analysis of Deviance Table

Model: binomial, link: logit

Response: Estado

Terms added sequentially (first to last)

        Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                      358     476.38              
Long     1  103.010       357     373.38 < 2.2e-16 ***
Sp       2   31.771       355     341.60 1.262e-07 ***
Long:Sp  2    0.459       353     341.15    0.7949    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interacción no significativa

glm_algas2 <- glm(Estado ~ Long+Sp, algas, family = "binomial")
anova(glm_algas2)
Analysis of Deviance Table

Model: binomial, link: logit

Response: Estado

Terms added sequentially (first to last)

     Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                   358     476.38              
Long  1  103.010       357     373.38 < 2.2e-16 ***
Sp    2   31.771       355     341.60 1.262e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(D2 <- 1 - (341.6/476.38))
[1] 0.2829254

Supuestos del modelo

par(mfrow=c(2,2))
plot(glm_algas2)

Estructura del modelo

summary(glm_algas2)

Call:
glm(formula = Estado ~ Long + Sp, family = "binomial", data = algas)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.80703    0.42701  -6.574 4.91e-11 ***
Long         0.26505    0.03174   8.351  < 2e-16 ***
SpTom       -0.57476    0.34528  -1.665    0.096 .  
SpVer       -1.84877    0.36933  -5.006 5.57e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 476.38  on 358  degrees of freedom
Residual deviance: 341.60  on 355  degrees of freedom
AIC: 349.6

Number of Fisher Scoring iterations: 5

Gráfico de predicciones

x_seq <- seq(min(algas$Long), max(algas$Long), length.out = 100)
newdata <- expand.grid(Long = x_seq, Sp = levels(algas$Sp))
# Obtener predicciones en escala de respuesta (probabilidades)
pred <- predict(glm_algas2, newdata, type = "link", se.fit = TRUE)
# Calcular intervalos de confianza
crit <- qnorm(0.975)  # 95%
newdata$fit_link <- pred$fit
newdata$lwr_link <- pred$fit - crit * pred$se.fit
newdata$upr_link <- pred$fit + crit * pred$se.fit
# Convertir a escala de probabilidad (logística)
inv_logit <- function(x) exp(x) / (1 + exp(x))
newdata$fit <- inv_logit(newdata$fit_link)
newdata$lwr <- inv_logit(newdata$lwr_link)
newdata$upr <- inv_logit(newdata$upr_link)
# Promediar sobre el factor 
pred_marginal <- aggregate(cbind(fit, lwr, upr) ~ Long, data = newdata, mean)
# Graficar relación promedio con IC
library(ggplot2)
ggplot(pred_marginal, aes(x = Long, y = fit)) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), fill = "#0072B2", alpha = 0.2) +
  geom_line(color = "#0072B2", linewidth = 1.3) +
  labs(
    x = "Talla individuo (cm)",
    y = "Predicción probabilidad estructuras reproductoras",
  ) +
  theme_minimal(base_size = 14)
newdat <- expand.grid(
  Long = seq(min(algas$Long), max(algas$Long), length = 100),
  Sp = levels(algas$Sp)
)

# Predicciones de probabilidad
pred <- predict(glm_algas2, newdat, type = "link", se.fit = TRUE)

# Convertir de escala logit a probabilidad
newdat$fit <- plogis(pred$fit)
newdat$lwr <- plogis(pred$fit - 1.96 * pred$se.fit)
newdat$upr <- plogis(pred$fit + 1.96 * pred$se.fit)

ggplot(newdat, aes(x = Long, y = fit, color = Sp, fill = Sp)) +
  geom_line(linewidth = 1.2) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.15, color = NA) +
  labs(
    x = "Talla individuo (cm)",
    y = "Predicción probabilidad estructuras reproductoras",
  ) +
  theme_minimal(base_size = 13)

Comparaciones entre especies

library(multcomp)
summary(glht(glm_algas2, linfct = mcp(Sp = "Tukey")))

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: glm(formula = Estado ~ Long + Sp, family = "binomial", data = algas)

Linear Hypotheses:
               Estimate Std. Error z value Pr(>|z|)    
Tom - Fra == 0  -0.5748     0.3453  -1.665 0.217830    
Ver - Fra == 0  -1.8488     0.3693  -5.006  < 1e-04 ***
Ver - Tom == 0  -1.2740     0.3126  -4.076 0.000131 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)